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Abstract 

In this manuscript, eigenvalues of the Electron Temperature Gradient (ETG) modes and Ion 
Temperature Gradient (ITG) modes are determined numerically using Hermite and Sine differenti- 
ation matrices based methods. It is shown that these methods are very useful for the computation 
of growth rates and threshold of the ETG and ITG modes. The total number of accurately com- 
puted eigenvalues for the modes have also been computed. The ideas developed here are also of 
relevance to other modes that use Ballooning formalism. 
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I. INTRODUCTION 



Computation of growth rates of microinstabilities in tokamaks is of considerable inter- 
est to the controlled fusion community. Often nonlocal analysis of Braginskii, Gyrokinetic 
or Gyrofluid equations leads us to a system of differential equations or differentio-integral 
equations with an unknown 'eigenvalue'. Such eigenvalue problems are encountered exten- 
sively in areas of ballooning modes, drift waves and other microinstabilities. Such equations 
have usually been solved by the shooting method, integral eigenvalue codes, phase integral 
methods and so on by a number of authors and Refs is certainly an incomplete list. 

However shooting methods suffer from a number of disadvantages. Shooting methods are 
very sensitive to initial conditions. It may be difficult to arrive at the right answer unless a 
very good approximation to the correct result is known. Furthermore, if determination the 
full eigenmode spectrum is intended then the whole 2-D parameter space of initial guesses 
Ur + must be spanned. Keeping in mind sensitivity to the initial guess, this may not be 
easy. Shooting methods can run into difficulties if the eigenvalues of the Jacobian matrix 
are widely separated. Often results are determined by round off errors rather than by the 
equation itself. Integral eigenvalue codes like those that use methods of Delves and Lyness 
or Davies Method also require initial guesses [2T] . 

Two methods that overcome many of these disadvantages are the Finite Difference (FD) 
and Spectral collocation methods Finite difference methods consist of replacing the 

derivatives by an approximation ( such as du/dx ~ [u{x+h)—u{x — hy\/2h where h is the grid 
size. This introduces truncation error and accuracy may be compromised (although a finite 
difference method is easy to program). Higher order finite difference methods impose strong 
stability restriction. Gary and Hegalson [19j have highlighted the usefulness of difference 
methods in computing eigenvalues of differential equations. Specifically, it has been pointed 
out by them that the shooting methods that use MuUer method or Lagurre method may 
run into difficulties (see table VI of Ref [12]). 

Spectral methods involve expansion of the solution by a finite sum in a orthogonal basis 
functions and minimizing of the residual [17] . Different spectral and pseudospectral methods 
differ in their minimization strategies. In this manuscript, the use of Differentiation Matrix 
(DM) based spectral methods to compute the ETG growth rate has been illustrated. Dif- 
ferentiation matrix based methods score over finite difference methods in their ease of use. 
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greater accuracy and faster execution for the same matrix size. The importance of these 
methods was reahzed by the author through an informal comparison observed in course of 
this work: a 15 x 15 DM gave similar accuracy as a 150 x 150 FD matrix. Since compu- 
tational time scales very fast (It costs ~ lOA^^ operations for the QZ [22l [23] algorithm), 
the ratio of the time consumed in the FD method (A^ = 150) while that in DM (A^ = 15) 
would be (150)^/15'^ = 1000!! - a huge increase in speed!. It must be noted however, the 
matrices involved in the DM and FD methods are, respectively, banded and full. This might 
make the DM method much more expensive for the same N number. Moreover, the FD 
performance could probably also be improved by optimizing the non-equidistant grid. 

Furthermore, as will be shown in this report, ETG growth rates converge faster with 
differentiation matrix methods. Differentiation matrix methods have one more advantage. 
The process of differentiation can be done through the matrix- vector product = D^^"* f 
where / is a vector of function values at the nodes {x^}. For example consider the model 
problem 

d^f/dx^ - df/dx = Xf 

which in matrix form becomes 

where D^^^ and Z^^^^ are a first and second order differentiation matrices respectively. Thus 
in these methods, the process of differentiation can be carried out through a matrix-vector 
product. Thus in principle one could create differentiation matrices for the curl or a gradient 
operator to solve more complex problems. The application of differentiation matrices is 
relavant to plasma physics and specifically to toroidal version of drift instabilities employing 
the ballooning formalism. 

It must be pointed out that the goal of this manuscript is not to show that finite differ- 
ence algorithms are not useful. In fact, a variety of situations can be envisaged where finite 
difference algorithms may be indispensible. Furthermore, neither the finite difference algo- 
rithm nor the differentiation matrix based methods used here are the very best and therefore 
a full academic comparison therefore cannot be done. Since a general analytical solution of 
the equations discussed here in not available, a comparision of numerical solutions from FD 
and DM methods is used as a means to fill this gap. 

The main goal therefore, is to present application of differentiation matrices to linear 
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analysis of ballooning modes. In the present paper, features of the collisionless, linear 
electromagnetic ETG mode[Tni|T2] are investigated used as a specific expample of ballooning 
type modes. The ETG mode is widely recognized as a mechanism for the small scale 
turbulence driven electron heat transport in tokamaks. The ETG mode is a short wavelength 
{pi » A_L >> pe), low frequency {u « Qe) fast growing mode {jetg - ^/e~r]^Ce/ Ln) 
driven unstable by temperature and density gradients. Here, pi is the ion gyroradius, pe is 
the electron gyroradius, Qe is the electron cyclotron frequency, e„ = 2L„/i?, rj^ = Ln/Lxe 
where L„ is the background density scale length. 

The organization of this paper is as follows. In the next section (Sec|ll]), the basic equa- 
tions of the ETG mode will be presented. An approximate solution in the strong ballooning 
limit is also derived and approximate eigenvectors computed (Figure [T]). Finite difference 



and spectral collocation methods form the subject of the sections III and |IV[ Section III 



presents the finite difference method used for benchmarking. An optimised finite difference 
method is also presented. Next, in Sec. Ill, comparisons with the 'analytical approximation' 
(Figure [T]) is also presented. Disadvantages of the method are also highlighted. In section 



IV, the differentiation matrix based method are compared with the finite difference method. 



For this purpose, comparison is made in three different ways: (i) Growth, (ii) Matrix size, 
and (iii) Conjugacy of eigenvalues. Next, the problem of determination of the number of 
eigenvalues is addressed for which the 'nearest' and 'ordinal' distances have been computed. 
A summary and conclusions are presented in Section |V| 



II. THE BASIC EQUATIONS OF ELECTRON TEMPERATURE GRADIENT 
MODE 

In this section the fluid description of the ETG mode will be discussed and a set of 
basic equations will be presented. An advanced electron fluid model for the core plasma (i.e. 
collisionless plasma) has been used. This fluid model is analogous to the fluid equations used 
in Refs. ([9], [10], [H])- The governing model equations are electron continuity equation, 
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(1) 



parallel momentum-equation and temperature equation: 

a; [r - Vi] + [1 - + r) + (1 + r],)Vl] ke4> 

i-il - Vi)zu; + f (1 + Ve)^ke)A^^ + ^(0 - n - f ) = 0, 
-uf + kelsnf + [(r/e - 2/3) A;e - f u;] = 0. 

Where fi, 0, T and A\\ are the normalized perturbations in density potential 
(50) , electron temperature (5Te) and parallel vector potential {eA\\/Te) respectively: 
[(j),h,f,6B\^ = L„/pe (e(50/Teo,5n/no,5Te/Teo,(55|j/5), A\\ = {2ceLn/ l3eCPe)eA\\/TeQ. The 
perpendicular scales are normalized to pe, the electron larmor radius while the par- 
allel scale is normalized to L„. As in the standard ballooning formalism, V^/ = 
-klf = -klf (l + {s9 - asmOf^ = e„ [cos^ + (s^ - asin^) sin^], Vy/ = ik\\f ~ 
(l/qR) {df/d6) where 6 is the extended coordinate in the ballooning formalism. Here u is 
the unknown, complex eigenvalue to be determined and e„ = 2Ln/R, where R is the major 
radius and Ln is the density scale length. Also note that this is a set of coupled, complex, 
differential and algebraic equations. 

First, an approximate analytical solution for a qualitative comparison is presented. A 
general analytical solution of Eq. [T]is not available, however the most important case where 
one can obtain an analytical solution is the strong ballooning limit. In the strong ballooning 
limit, one can assume that g{6) = cos{6) + {s9 — asin{9)Y is slowly varying with 9 so that 
g{d) ~ 5((0) = 1. In this limit, Eq. [l| can be written as 
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where Tq = kjXj^^ + t + 6, u = coLn/ce and k = kop^. One specially looks the solutions of 
Eq. [2] of the form 

$(^) = exp(-(7^V2), (6) 



where Re{a) > 0, where a = ±^/C/A, this gives us the dispersion relation 



B = -a A. (7) 



Here A, B and C are parameters independent of d as in Ref [TT]. Note that the solutions of 
the form $(^) oc F{6)exp{~cr6'^) were used. Especially the case F{6) = 1, which corresponds 
to the above solution. 

Next, a discussion of the numerical solution of the equations [sjjT] will be discussed. This 
will be used as a verification benchmark for the finite difference and differentiation matrix 
based methods of solution of Eqs. ([T]). The function h(uj) = B -\- a A is a nonlinear function 
of the complex variable oj = cor + icoi, hence we have solved Eq. ([t]) by finding the roots of the 
function h{uj) using Muller method. In figure [l| the eigenvectors corresponding to F(9) = 1 
from the above 'analytical' solution have been drawn as solid lines. The dots in this figure 
will form the subject of our next discussion - the finite difference method, which is discussed 



Sec. Ill The solid squares in this figure refer to computation using the Hermite Method as 



descibed in Sec. IV, Please note that in this figure, for easy comparison the eigenvectors 
have been normalized in such way so that Real[^{6 = 0)] = 1 and Imag[^{d = 0)] = 1 for 
all the methods. 



III. FINITE DIFFERENCE METHODS 

The basic equations solved in this manuscript, Eq. ([T]), can be considered as a system of 
first order ordinary differential equations of the form 

^ + J] 6, {x)yj = A J] a,- (x)?/,- (8) 

Often, in the past, a number of authors have solved similar eigenvalue problems but as 
a single second order differential equation that involves powers of the eigenvalue A in the 
following manner: 

A) ^ = ^ b{x, \)y (9) 

where a(x. A), A) can be nonlinear functions in x but polynomial in the eigenvalue A. 

Finite difference methods are the most widely used methods to solve for eigenvalues of 
such ODE-eigenvalue problems. The basic ideology that we have used is highlighted in the 



6 



following equations: 



Vk - Vk-i 




(10) 



An implementation of such finite difference methods is explained in Ref [5^. Although 
the method is elaborate, it has the following disadvantages. This method computes the 
eigenvalues by calculating the zeros of a determinant. These zeros are computed using an 
iterative process which in turn, requires guesses. Requirement of guess values that makes it 
no better than shooting methods discussed earlier. Furthermore, for efficient use, the basins 
of attraction of the various initial guesses must be known beforehand. 

Instead of using the above code, a finite difference code has been developed to overcome 
the above shortcomings. This code has the following advantages not present in the above 
code: (i) generates the full eigenvalue spectrum simultaneously; (ii) there is an option to 
refine a particular eigenvalue; and (ii) for small changes in a parameter, the changes in 
a branch can be tracked automatically, otherwise it requires careful human intervention. 
Other salient feature are: (i) Can solve complex coupled differentio-algebraic equations; and 
(ii) arbitrary mesh (which is necessary and useful since boundary conditions are at infinity 
and the eigenfunction is expected to be localized). 

Next, the finite difference numerical scheme used is explained in detail. First, an odd 
number of grid points such that = is the central point lo is chosen. Next, allocate grid 
points such that 6{Io + j) = —0{Io — ])■ The points 6{I) need not be equally spaced but can 
be decided according to the problem at hand. The allocation of grid points was not optimized 
but chosen intuitively such that they are either (i) equally spaced; or (ii) 9 oc F{9i) where 

6i = 0, 69, 2 69, , 0.99 is a set of equally spaced points between and 0.99. The function 

F{9i) can be any monotonically increasing function of 9i like 1/(1 — ^i), exp{c9i) and so 
on. Next one discretizes the two differential equations at half grid points and the algebraic 
equation at full grid points. All eigenvalues and eigenvectors can be obtained simultaneously 
if the resulting set of coupled algebraic equations are solved by the QZ algorithm. 

In the strong ballooning limit, the full numerical computation closely follows the analyt- 
ical results. In Figure [T| eigenvectors from both the analytical and numerical computations 
are depicted as a validation of the code. The solid lines indicate analytical result while the 
solid circles depict the numerical values. It must also be pointed out that in this figure, the 
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X-axis range has been reduced (from the actual range used in the computation) to show that 
the eigenvectors from the analytical approximation closely follow the numerical result. The 
non-uniform distribution of grid points just explained is also illustrated in Figure [T] as solid 
dots. 

This finite difference method is not without disadvantages. Eq. ([T]) has three equations 
(two ODE's and one algebraic equation). Thus even for a modest mesh size of 150 points 
generates two matrices each of size (3 * 150) x (3 * 150). Investigation of the ETG mode 
requires variation of the basic parameters over a wide range. For example, for each value of 
shear parameter s and ballooning parameter a, safety factor q, plasma beta /3, temperature 
ratio r, eigenvalues for at least 200 — 300 values of rji must be computed. This large number 
of values of rji are used to accurately compute the threshold and this increases the time 
taken by the computer. 



IV. EIGENVALUES USING DIFFERENTIATION MATRICES 

Differentiation matrices are novel, powerful methods of finding the eigenvalues of differ- 
ential equations. These methods are very accurate. For smooth problems, convergence rates 
of the order of 0{e~^^) are often achieved [71 [8] as compared to finite difference methods 
where the convergence rate are slower, typically (A^~^) or (A^~^). 

The method employed is briefly explained below: A spectral collocation method fTU] for 
solving differential equations that is based on weighted interpolants of the form: 

has been used. The points {xj}^^^ is a set of distinct interpolation nodes, a{x) is a weight 
function, fj = f{xj) and the set of interpolating functions {(pj{x)}jL^ satisfies (pj{xk) = Sjk- 
Where 6jk is the Kronecker delta function. The differentiation matrix is then defined as a 
matrix with the entries: 

The process of differentiation can be done through the matrix- vector product = D^^^ f 
where / is a vector of function values at the nodes {xk}. 

This manuscript, mostly uses the Hermite Collocation method. In the Hermite Collo- 
cation method, the nodes Xi,X2,X3, ....x^ are the roots of H]^{x), the Hermite polynomial. 



The weight function is the Gaussian function a{x) = exp{—x'^/2) and the interpolant is 
defined as 



= . . (14) 



where, 

H'^{xj){x - Xj)' 

It must be noted that the real hne (— oo, cxd) may be mapped to itself by a change of variable 
X = bx, where b is any positive real number. Thus the freedom offered by the parameter 
b may be exploited to optimise the accuracy of the Hermite differencing process. Before 
proceeding further let us discuss the solution of Hermite differential equation itself. The 
solution the Hermite differential equation using a Hermite DM was performed. The errors 
in the first twenty five eigenvalues was found to be less than 10~®. 

The application of this method to plasma physics is not new. It has been used in sim- 
ulations of Vlasov equation since the sixties. Low order Hermite series (A^ = 2, 3) were 
shown to suppress numerical instabilities in a Vlasov simulation by Engelmann et al. |20] . 
More recently Shumer et al. have used velocity-scaled Hermite representations [18]. In this 
manuscript, we employ Hermite differentiation matrices to compute eigenvalues of the ETG 
instability. 



A. Comparison with FD: Growth Rate 

The eigenvalue of the ETG mode has also been computed using the following basic 
parameters s = 0.8, r = l,kg = 0.6, q = 1.4:, en = 0.509 and (3 = 0.01, while rji is varied. A 
typical result from these simulations is depicted in Figure |2j In this figure, the results of 
computation from Finite difference method is shown as solid lines and that from Hermite 
method is shown as circles. As can be seen in the figure, the Hermite method requires only 

= 16 to achieve accuracy comparable with a Finite difference computation with N = 121. 
The finite difference code uses an odd number of points 51, 121 etc. so as to have equal 
number of points to the left and to the right of ^ = 0. Since the Hermite method requires 
only a small number of points, its is faster by orders of magnitudes. This allows greater 
parameter scans. 
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B. Comparison with FD: Matrix Size 



The eigenvalues of the ETG mode for different values on N using the finite difference 
method and Hermite spectral collocation method have also been computed. As depicted in 
Figure ([s]), the Hermite method converges much before the two finite difference methods. 
In Figure ([s]), two finite difference methods are illustated by solid triangles and crosses. 
These FD methods only differ in terms of how the mesh was optimized. Below = 20, 
the finite difference code (which is illustated by triangles) gives wrong results {i.e., ~ 0), 
slowly this FD result begins to converge and a stable value is obtained only beyond = 100. 
To improve the convergence of the FD method, the mesh was optimized so that accuracy 
improves for A^ < 20. This is shown as crosses in Figure ([s]). The Hermite method, however 
gives us an acceptable eigenvalue even for A^ = 10, or A^ = 15. Note that the computational 
time taken for an accurate computation would also be small since A^ is small. 

The rate of convergence of these three methods is illustated in Figure Q where the error 
S{N) is plotted versus A^. Here, 

u.iN)-u^{N = 100) 
= (AT = 100) 

and u:^{N = 100) is the growth rate from the Hermite method at A^ = 100. In this figure, 
6{N) from the unoptimised finite difference method is represented using solid triangles, S{N) 
from the optimised finite difference method is represented by crosses and the S{N) from the 
differentiation matrix (Hermite) method is represented by solid dots. The solid line in this 
figure represents an error level of 1% which is easily achieved by the Hermite method. The 
error also decreases fastest for the this method. 



C. Comparison with FD: Conjugacy of Growth Rates 

Looking at the basic equations, one realizes that for a particular set of parameters, the 
maximum growth rate {'jmax) also has a complex conjugate 7mm- Mathematically, one 
expects that •jmax = —7mm- This however may not be true in reality. Only one of the two 
roots may be more accurate than the other. Thus the difference between the two growth rates 
can be considered to be a measure of the numerical errors in computation. The parameter 
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(^7 where, 

57 = , (16) 

'~1max '~1min 

can be considered as a measure of this error. To illustrate the competitiveness of the method, 
the variation of ^7 for various N , using the finite difference as well as the Hermite Method 
has also been computed. This is shown in Figure [5] where the triangles refer to a Hermite 
calculation while the dots refer to finite difference method. Values from both methods are 
small but the Hermite method appears to be slightly more accurate than finite difference. 
It is possible that the slightly smaller inaccuracies of the Hermite-case are simply due to a 
slightly unfavorable ordering in the matrix or round off errors. It must also be noted that 
this error is also smaller compared to the error in the approximation schemes. 



D. Number of Convergent Eigenvalues 

The system of equations [T] has three equations in three variables. Increasing N, further- 
more, would give us 3A^ eigenvalues. For an arbitrary N, increasing to infinity, the number of 
eigenvalues would also be large. All of these eigenvalues might not be good approximations 
to the eigenvalues of the actual differential equation. Physically one expects only a few of 
these should remain the same for every N. Thus it is important to ascertain the number of 
true eigenvalues of the differential equation at hand. 

The following two situations can be envisaged: (i) the eigenvalue 'drifts' as N increases, 
(ii) a new (probably spurious) eigenvalue is inserted between the two earlier ones. In order 
to ascertain the correct eigenvalues, Boyd et. al. have proposed a scheme [6] for ocean tides 
bounded by meridians. Following this prescription, the 'ordinal' and 'nearest neighbor' drifts 
can be computed as explained below. 

One computes eigenvalues using two different matrix sizes A^i and N2 such that N1/N2 > 
1/2. The physically important eigenvalues include the highest growing {'~fmax > 0) and the 
most stable (7mm < 0) These modes are often complex conjugates of each other {jmax ~ 
—Imin) ■ Thus growth rate spectrum is then sorted in two different ways- ascending and 
descending. 

The ordinal drift is defined by 

Sj,ordinal = ~ ^f'^l^ (!''') 
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where A^^ is the imaginary part of the jth eigenvalue (after the eigenvalues have been sorted) 
as computed using N spectral coefficients. 
While the nearest neighbor drift is given by 

next, the intermodal drift is defined as follows 



Af-Afi, zfj 



HlAj - Aj-il + - Xj\ if j > 1, 
Boyd et. al. then recommend making a log/linear plot of the following ratios: 

_ mm(|Aj|,a,) 

' j, ordinal i- \ 

('j, ordinal 

_ min{\Xj\,aj) 

' j, nearest r- \^^/ 

^j,nearest 

which must both be large for the eigenvalue computation to be accurate. Boyd et. al, 
illustrate for Legendre's differential equation in Ref [B], Figure 4. 

The ordinal and nearest distances for the ETG mode has been computed using the pa- 
rameters s = 0.8, r = l,kg = 0.6, g = 1.4, e„ = 0.909, /3 = 0.01 in Equation [l]. The choice 
of parameters is not arbitrary, they have been chosen to be the same as those of Jenko et. 
al. [21]. Next, the eigenvalues are sorted in two different ways- ascending and descending 
orders. A total of ten numerically accurate eigenvalues were found, five each for each case 
signifying five growing modes and five damped modes. The case of highest growing modes 
is illustrated in figure[6]. Such a large number of eigenvalues can be understood using the 
following argument. The argument is based on the results of the local and nonlocal analy- 
sis of the electron temperature gradient modes given in Ref. [H] by R. Singh et. al. The 
Equation 31 of Ref [11] is a third order dispersion relation retaining parallel dynamics. That 
is one growing, one stable and and real frequency branch. Moving to nonlocal. Equations 
41-44, we realize that the dispersion relation is now a fifth order. Furthermore Eqns 41-44 
assume solutions of the type exp{—a6^/2) and ignores the higher branches with the form 
Hn{6)exp{—cr6'^/2), n = 1,2,3 where Hn is the Hermite polynomial. Higher values of n 
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add more branches to be found numerically. Thus the possibility a large number of eigen- 
values cannot be ruled out. This conclusion is not new. Gao et. al. [13j have identified 
a series of unstable branches for the ETG mode. Similarly (THj et. al. have shown six 
growing modes. Lee et. al [[25j ]have studied the three most unstable modes as a function 
of parameters. 

One might also ask oneself if the number of these eigenvalues changes with physical 
parameters. The variation of the number of good, highest growing eigenvalues with the 
parameter e„ has aslo been investigated. For large values of e„ (~ 0.9), the number of 
good eigenvalues is five. However on decreasing e„, (~ 0.5), the number of good eigenvalues 
increases and becomes seven. This is illustrated in figure [7j On a further decrease of e„, 
(~ 0.2), the large separation between good eigenvalues and spurious diminishes with the 
mode number j. 



E. The Sine Method 



The use of Sine method to compute the eigenvalues of the ETG mode has also been in- 
vestigated. The Sine Method is usually applicable on an interval (— oo, oo) using equidistant 
nodes with spacing h and the nodes 

+ 1 

Xk = {k —)K k = l,...,N (21) 

and the weight functionQ;(x) = 1 and the interpolant 

Like the Hermite method discussed earlier, the Sine method also contains a free parameter 
h. The results are encouraging and represented in Figure [8j It gives a reasonably accurate 
growth rate for even N=3. Thus the Sine method is very useful if one is interested only in 
the computation of growth rates. The real frequency for = 3 is wrong by about 20% but 
rapidly improves as increases. 



F. Computation of the ITG growth rate 

Apart from investigating the ETG mode in detail, some computations of the ITG mode 
using a different set of basic equations has also been attempted. The following are not a set 
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of differentio-algebraic equations as Eq. ([T]), but a single second order differential equation 
as derived in Ref. [H]: 

^5 = h[in- 1)A + 2eng{e) + hn] ij (23) 
du 

where 

Where ip is the eigenvector proportional to the normalized eigenvector e(f)/Te as in Ref 
[H] . The parameter 6 is the extended ballooning coordinate, q is the safety factor and 
s = {r/q)dq/dr is the shear parameter, g{6) = cos{6) + s6sin{6) is a function which will 
be symbolized by g from now on. The above equations can be converted into a polynomial 
eigenvalue problem of the form "^f^QAp^fl^X = or 

+ [C5(-Ci + CsOij + C2^]n (25) 

+[C5(-i + Ci + + C3b)^p]n^ + [C5(i + bC2)^] n^ = o 

where 6*1,6*2. ..C5 are independent of fl and are defined as follows: 



Ci = ^ C2 = l + ^ (26) 



T 6T 

n 

The results of the numerical computation are depicted in Figure |9] where the dashed line 
corresponds to the growth rate and the dashed dotted line corresponds to the real frequency. 
The parameters used are e„ = 1, r^j = 8, g = 1, r = 1, 6e = 0.1 and s = 1 and the boundary 
condition ip = at x = ±00. As in the ETG case, the operator d'^ip/dO'^ is replaced by the 
corresponding Hermite differentiation matrix. 



V. SUMMARY AND REMARKS 



Two methods to numerically solve for eigenvalues of the Electron Temperature Gradient 
(ETG) mode were presented. These methods were Hermite and Sine methods. For verifi- 
cation purpose the growth rate was also computed analytically as well as numerically using 
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finite difference metfiods. Tlie advantages and disadvantages of eacli of tliese metliods were 
also liigliliglited. 

Altliougli tliese methods seem promising they also suffer from some disadvantages. Firstly 
banded matrices of finite difference matrices are replaced by full matrices. Secondly the 
scaling parameter b in the Hermite method or the spacing h in the Sine method must be 
adjusted for optimum results. 

For computing ETG eigenvalues, spectral collocation methods were found to be more 
accurate, faster and more useful than finite difference methods. As another example, the 
eigenvalue equation for the Ion Temperature Gradient (ITG) modes was also attempted. 
Furthermore, recenly, these methods were also applied to the Drift Resistive Ballooning 
Modes (DRBM) with encouraging resultspGj. Finally, the results presented here are quite 
general and can be applied to a number of other microinstabilities. 
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FIGURE CAPTIONS 



FIG 1. Comparison of eigenvectors from different methods is shown. The analytical 
eigenvector from Eq. (|6|) is shown as a solid line and the numerical computation using 



the finite difference method as described in Sec III is shown as dots (and labeled FD in 
legend). The solid squares refer to computation using the Hermite Method as descibed 
in Sec. IVl 

FIG 2. Real frequency Ur and growth rate Ui from finite difference method is compared 
with Ur and Ui from Hermite DM method. Note that Hermite Method with N = 16 
compares well with the finite difference method at = 121 when Ui > . 

FIG 3. Real frequency Ur and growth rate Ui from the optimsed finite difference 
method (indicated by crosses) is comapred with the finite difference method (triangles) 
and the Hermite DM which is indicated by solid dots. 

FIG 4. Numerical convergence of ETC growth rate as a function of N. Here, the error 
6{N) = [cOi{N) -co^iN = 100)]/u^{N = 100) is plotted versus N for three diflFrent 
methods finite difference (triangles), the optimsed finite difference method (crosses) 
and the Hermite DM method (dots). 

FIG 5. Numerical convergence of ETC growth rate is illustated through the congugacy 
of growth rates with the help of the ratio ^7 . 

FIG 6. Ordinal and nearest distances plotted on a logarithmic scale versus mode 
number j for the ETC mode, e„ = 0.9. 

FIG 7. Ordinal and nearest distances plotted on a logarithmic scale versus mode 
number j for the ETG mode, = 0.5. 

FIG 8. Solution of the equations for the ETG real frequency Ur and growth rate Ui 
using Sine Method is compared with the finite difference method. Note the small value 
of = 3 compares well with A^ = 121. 

FIG 9. Growth rate (solid line) and real frequency (dashed-dotted line) of the ITG 
mode as a function of e„. 
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